
function [conv_crit, A, K] = tremble_r_fp(r, param, util_type, N, smooth_wind, T, yi, reset_shocks, stateGrid_size,eps_eta_mat,sigma_tremble, p_unemp, w_unemp)

bta = param(1);
abar = param(3);
W_cc = param(4);
kappa = param(5);
deltax = param(6);
sigma_c = param(7);
psi = param(8);
cov_fun = param(9);
gam = param(10);
sigma_l = param(11);
lbar = param(12);
delta_eta = param(13);
delta_inh= param(14);
delta    = param(15);
alpha    = param(16);
b        = param(17);
btat     = param(18);
amin     = param(19);
amax     = param(20);
na       = param(21);
crit_dist = param(22); 


% Step 1: get aggregate capital and wage implied by r
K = ((r + delta)/alpha)^(1/(alpha-1)); % Aggregate demand of capital (from firm k-FOC)
w = (1-alpha).*K.^alpha; % Wage (from firm l-FOC)

% Step 2: get optimal policy function

rt = (r*(1-deltax*delta_inh) + deltax*(1-delta_inh))./(1-deltax); % effective interest rate (includes insurance market)

[xs,apstar,cstar,~,~] = bc_opt_pol_iid_egm(btat, rt, b, log(w) + lbar, sigma_l, gam, amin, amax, na, util_type, p_unemp, w_unemp);
[xs,xi] = unique(xs);
temp = griddedInterpolant(xs,cstar(xi),'linear','linear');

%stateGrid_size = 1001;
stateGrid = linspace( xs(1), xs(end),stateGrid_size);
stateGrid_step = stateGrid(2) - stateGrid(1);
cstar_fixedpoints = temp(stateGrid);

temp = griddedInterpolant(stateGrid,cstar_fixedpoints,'linear','linear');
cstar = @(x) temp(x);
%apstar = @(x) fast_interp( x, apstar_fixedpoints, stateGrid_step, stateGrid,stateGrid_size);

% Step 3: compute stationary distribution of a
ap = 0.*ones(1,N); % initial distr of wealth

%conv_distr = 0;

conv_mat1 = zeros(1,smooth_wind); % keeps track of the mean for each of the last 30 iter
conv_mat2 = ones(1,smooth_wind); % keeps track of the mean for each of t=end-30 to end-60


for t = 1:T
    
    y = w*yi(t,:);
    
    a = ap.*(reset_shocks(t,:)>deltax*(1-delta_inh));
    
    ap = (y + (1+rt).*a - max( [cstar(y + (1+rt).*a) + eps_eta_mat(t,:)*sigma_tremble;1e-3*ones(size(y))] )); % a is multiplied by 1 - death indicator
    
    conv_mat_index = mod(t,smooth_wind)+1;
    conv_mat2(conv_mat_index) = conv_mat1(conv_mat_index);
    conv_mat1(conv_mat_index) = mean(ap);
    
    %conv_distr= abs(mean(conv_mat1)-mean(conv_mat2));
    
end

% Step 4: compute aggregate assets by integrating the distribution of a

if  abs(mean(conv_mat1) - mean(conv_mat2)) > crit_dist
    disp(['r = ', sprintf('%1.5f' , r), '     !Distribution did not converge!']) 
else
    disp(['r = ', sprintf('%1.5f' , r)]) 
end

A = mean([conv_mat1 conv_mat2]);

conv_crit = A - K;

end
